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Abstract 

Global hybrid (electron fluid, kinetic ions) and fully kinetic simulations of the magnetosphere 
have been used to show surprising interconnection between shocks, turbulence and magnetic re- 
connection. In particular collisionless shocks with their reflected ions that can get upstream before 
retransmission can generate previously unforeseen phenomena in the post shocked flows: (i) forma- 
tion of reconnecting current sheets and magnetic islands with sizes up to tens of ion inertial length, 
(ii) Generation of large scale low frequency electromagnetic waves that are compressed and ampli- 
fied as they cross the shock. These “wavefronts” maintain their integrity for tens of ion cyclotron 
times but eventually disrupt and dissipate their energy, (iii) Rippling of the shock front, which 
can in turn lead to formation of fast collimated jets extending to hundreds of ion inertial lengths 
downstream of the shock. The jets, which have high dynamical pressure, “stir” the downstream 
region, creating large scale disturbances such as vortices, sunward flows, and can trigger flux ropes 
along the magnetopause. This phenomenology closes the loop between shocks, turbulence and 
magnetic reconnection in ways previously unrealized. These interconnections appear generic for 
the collisionless plasmas typical of space, and are expected even at planar shocks, although they 
will also occur at curved shocks as occur at planets or around ejecta. 

PACS numbers: 52.35.Vd, 52.65.-y 
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I. INTRODUCTION 


There is growing evidence that magnetic reconnection and turbulence are intimately 
connected in magnetized plasmas (e.g., [1-5] and references therein). However, the link 
between shocks and these two processes has been less clear. We explore this connection 
using global kinetic simulations of the Earth’s bow shock, and in the process address a 
number of puzzling observations. Due to the availability of in situ measurements [6], the 
Earth’s bow shock is the most studied example of collisionless shocks. The bow shock is a 
curved shock formed due to interaction of the supermagnetosonic solar wind with the Earth’s 
magnetic dipole which acts as a barrier. The shocked solar wind is called the magnetosheath 
and is bounded by the magnetopause which diverts the flow. 

Observations, along with theory and simulations, have provided considerable insight into 
the physics of collisionless fast shocks. The curved nature of the bow shock means that 
at a given instant, there are regions with quasi-parallel Qy (0° < 9bn < 45° and quasi- 
perpendicular Q± (45° < 0 bn < 90°) geometries. Here 6bn is the angle between the shock 
normal and the incident magnetic field. In the quasi-parallel case, there exists a significant 
number of ions reflected from the shock that reach to large distances upstream. The relative 
streaming between these ions and the incoming plasma excites instabilities which generate 
low frequency waves that can grow to large amplitudes, giving rise to steepened fronts called 
magnetosonic shocklets and short large-amplitude magnetic structures (SLAMS) as well as 
other nonlinear structures such as cavitons (see review in [6] ) . These structures are convected 
back into the shock and contribute to generation of turbulence in the magnetosheath. 

There are also reflected ions associated with quasi-perpendicular shocks, but these ions 
remain confined to the foot of the shock, leading to shock reformation on the gyroscales of the 
ions for a certain range of parameters (e.g., [7]). There can also be turbulence downstream of 
quasi-perpendicular shocks due to temperature anisotropy driven instabilities with magnetic 
fluctuation levels B rms /B 0 ~ 0.1 as compared to B rms /B 0 ~ 1 downstream of quasi-parallel 
shocks [8]. 

Despite this progress, a number of important issues remain. One issue of considerable 
debate is whether bow shock turbulence has significant large scale consequences, especially 
since the turbulence is initiated on ion scales. In particular, observations have shown evi- 
dence of anomalous flows in the magnetosheath. Examples include observations of fast jets 
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(see [9, 10] and references therein) and sunward flows (see [11] and references therein). The 
origin of such flows and their impact on the magnetosphere remain controversial. Since solar 
wind has embedded discontinuities that frequently impinge on the magnetosphere, models to 
explain these anomalous flows fall into two general categories: those that invoke the role of 
discontinuities and those that rely on the physics of the bow shock itself. Similarly, observa- 
tions have reported reconnecting current sheets in the quasi-parallel magnetosheath[12, 13]. 
These are difficult measurements and confirmation through simulations would be helpful 
but no evidence for such current sheets have yet been reported in the global simulations. 

Simulations provide a global view of the magnetosphere and its dynamical evolution. In 
order to address the above issues, and to gain a better understanding of the turbulence 
associated with the bow shock, we use global hybrid and global fully kinetic simulations 
that capture the self-consistent formation of the bow shock and its associated turbulence. 
Most studies to date have focused on the upstream region of the Q\\ bow shock, called 
the ion foreshock. We use our simulations to compare and contrast the nature of fluctua- 
tions/turbulence in the ion foreshock, Q|| magnetosheath, and Q± magnetosheath and their 
large scale consequences. We have run the simulations to a factor of ~ 9 times longer 
than most previous studies in order to gain insight into the evolution of turbulence in time. 
Through these simulations, we find that the Q\\ magnetosheath turbulence is much stronger 
than either that in the foreshock or in the Q± magnetosheath. We demonstrate, for the 
first time, that this enhanced level of turbulence leads to generation of reconnecting current 
sheets and formation of magnetic islands. This finding establishes the link between shocks, 
turbulence and reconnection. This result is consistent with observations of reconnection in 
the quasi-parallel magnetosheath [12, 13], where, however, no search for magnetic islands 
was performed. In 3D, magnetic islands correspond to flux ropes. In the remainder of this 
paper, we will use the term magnetic islands and flux ropes interchangeably. Since some of 
the flux ropes reach tens of ion inertial length, we expect a dedicated observational search 
for flux ropes in the magnetosheath to yield many examples. One possibility to search for 
such flux ropes in the data may be through the use of the Grad-Shafranov reconstruction 
technique [14, 15]. 

The simulations also enable us to address large scale consequences of turbulence in the 
magnetosphere. We examined the possibility that foreshock turbulence, rather than inter- 
action of solar discontinuities with the magnetosphere, could cause anomalous flows. We 
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show that turbulence causes corrugation of the shock surface which in turn gives rise to 
collimated, high velocity plasma jets inside the magnetosheath. This is in agreement with 
the model originally proposed by Hietala et al. [9] and the properties of the jets are found 
to be consistent with high velocity magnetosheath flows recently reported in the observa- 
tions [9, 10]. The jets have high dynamical pressure and we show that upon reaching the 
magnetopause they can at times trigger so-called flux transfer events [16]. 

Next, we considered whether the jets, which penetrate into the magnetosheath, can trigger 
sunward flows. Such a connection was proposed by Shue et al. [11] who attributed the 
origin of the sunward flow to the interaction of a jet with the magnetopause, which in turn 
would cause the magnetopause to rebound, creating a sunward motion. However, we find 
sunward flows in the vicinity of many jets that terminated in the magnetosheath far from 
the magnetopause. A careful examination of the time history of the flow reveals a complex 
process that would have been difficult to uncover in observations. In this process, the jets 
“stir” the magnetosheath, creating large scale disturbances and pushing the ambient low 
density plasma and the embedded magnetic field. The plasma in the sunward flow regions is 
not the plasma in the jets that was deflected, rather it is the plasma accelerated by thermal 
and dynamical pressure that develop because of the jets. In other words, the energy of the 
jets is convected into sunward flow rather than the deflection of the jets themselves. This 
energy conversion can give rise to sunward flows that extend tens and even hundreds of 
ion inertial length into the solar wind, well past the nominal position of the bow shock. 
The strong turbulence in the magnetosheath also gives rise to regions where the magnetic 
field during anti-sunward IMF orientation can point sunward in the magnetosheath. This is 
another testable prediction for observations. 

Another consequence of the jets is found to be the nonlinear formation of large scale 
vortices in the magnetosheath. These vortices are driven due to turbulence and are to 
be distinguished from the usual Kelvin- Helmholtz driven vortices along the flanks of the 
magnetopause ([17, 18] and references therein). 

Although these results were obtained from simulations of the bow shock, a curved shock, 
similar effects are also expected at planar shocks. For example, the current sheets formed 
downstream of the quasi-parallel shocks is due to the strength of the turbulence, rather than 
any physics tied to the curvature of the shock. Similarly, the rippling of the surface of the 
shock is due to interaction of the shock with the steepened fronts getting convected back into 
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it, an effect that is also present in planar shocks. The paper concludes with the discussion 
of the implications of these findings for the general understanding of plasmas and draws the 
link between shocks, reconnection and turbulence. 


II. SIMULATION SETUP 

Our hybrid simulation code is H3D [19-21], a massively parallel hybrid code with a 
stretched mesh capability. Although we have conducted 3D global simulations of the mag- 
netosphere [20, 21], here we will focus on 2D simulations for three reasons: a) One of the 
key issues considered is whether the shock generated turbulence can lead to formation of 
reconnecting current sheets. The identification of reconnection in 3D remains a subject of 
active research and is still not completely resolved, b) Confirmation of the formation of 
reconnection sites observed in our global hybrid simulations through equivalent fully kinetic 
global simulations can only be done in 2D since 3D global fully kinetic simulations will re- 
main out of reach in the foreseeable future, c) We have conducted an extensive parameter 
study to identify the most interesting regimes which serve as conditions for 3D simulations 
to be published elsewhere. 

Distances are normalized to ion skin depth di, time is normalized to ion cyclotron fre- 
quency U ci , magnetic field is normalized to upstream magnetic field B 0 , density is normalized 
to upstream density N a and velocity is normalized to the Alfven speed Va = B a / (47T NoUii) 1 ^ 2 
where m ? ; is the ion mass. The ratio of plasma to ion cyclotron frequency is set as 6000 and 
the ions and electrons are taken to have the same temperature T, = T e for these runs. The 
Alfvenic Mach number is given by Ma = V sw /Va where V sw is the solar wind speed. We use 
a polytropic electron model with 7 e of 5/3 and 7/6, with the latter being suggested as more 
appropriate for shock physics[22, 23]. The deHoffman- Teller cross shock potential under 
polytropic electron closure is proportional to 7 e /(7 e — 1 )[kT e ]. Thus using the more nearly 
correct 7 e for solar wind thermal electrons [22, 23] of 7/6 leads to a larger deHoffmann- Teller 
potential in hybrid simulations, raising the potential coefficient to 7 versus 5/2 for y e of 5/3. 

We have conducted hybrid simulations for a range of D p of 40 to the Earth’s size of 600 
di, plasma beta of 1 to 4, Alfvenic Mach number Ma of 4 to 12, interplanetary magnetic 
field (IMF) angle of 0° to 90° and IMF B z /B a of 0 to 1. This range of parameters is not 
meant to be exhaustive but do present a good coverage of the relevant parameter space. 
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Here D p is the nose position of the magnetopause normalized to ion inertial length, which 
has been shown to be a good way to characterize the type of magnetospheric structure as a 
function of dipole strength [24], We have also examined the effects of i) spatial resolution 
by varying the grid from 0.25<i; to Id*, ii) particle per cell by considering a range of 10 to 
200, iii) and resistivity by considering different models such as no imposed resistivity, and 
several models based on gradients of magnetic field. We have found no significant effect of 
the resistivity on the physical processes considered here. However, care is required in general 
since use of large uniform resistivity as is done in some studies can lead to excess heating 
of the magnetosheath. We will focus this paper around the runs listed in Table 1 that are 
from converged solutions. 

The simulation setup is similar to that used in global MHD simulations and consists of a 
magnetic dipole embedded in a uniform solar wind plasma which is continually injected from 
the left boundary. Mirror dipole is used to ensure the magnetic field is just the IMF field at 
the injection boundary. This reduces magnetic field perturbations at the injection boundary. 
All other boundaries use open boundary conditions. As in global MHD simulations, an inner 
boundary is placed around the Earth to limit the Alfven speed from reaching relativistic 
values. Two particle boundary conditions are implemented: reflecting and absorbing. The 
dipole can also be ramped up in time, although we have not found significant differences 
compared to the case where the dipole is initialized at full strength at the start of the 
simulation. As the solar wind plasma interacts with the dipole, the magnetosphere forms in 
time. In 2D, the dipole is a line dipole in order to have a divergence free field. 


III. COMPARISON OF THE TURBULENCE REGIONS 

Figure 1 shows the intensity plot of magnetic field at three different times for run 1 where 
the IMF direction changes at the injection boundary from —45° to 10° at £l ci t = 90. This 
launches a rotational discontinuity (RD) that propagates towards the magnetosphere. The 
RD eventually runs into the reflected ions in the ion foreshock before reaching the bow 
shock. The time evolution of this interaction is shown in Figure 2 where we show a segment 
of the simulation zoomed in around the RD-reflected ion interaction region, along with ID 
cuts of several parameters at y — 3100^. Figure 2 shows evidence for particles following 
the field lines across the RD into the upstream region. An observational evidence for this 
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TABLE I: List of parameters for the global hybrid simulation runs shown in this paper. The 
number of particles per cell was 200. The spatial resolution was 0.5 dy except for runs 1, 5 (1 di) 
and run 6 (0.25 di). 


run 

%max 

Vniax 

1 

2048 

8192 

2 

2048 

4096 

3 

2048 

4096 

4 

8192 

8192 

5 

2048 

4096 

6 

512 

512 

7 

4096 

2048 


D p 

IMFB Z 

M a 

100 

0 

8 

100 

0 

8 

100 

0 

8 

300 

0 

10 

100 

0.6 

8 

50 

0.6 

8 

150 

0.6 

10 


IMFy 

imf 2 

tfiip 

-45° 

10° 

90 

10° 

10° 

N/A 

-45° 

10° 

90 

-45° 

10° 

150 

10° 

10° 

N/A 

10° 

10° 

N/A 

10° 

10° 

N/A 


effect would be particles on magnetic field lines that do not point at the bow shock. Figure 
2 also shows formation of at least three nonlinear structures: a fast shock, a cavity with hot 
plasma with depleted plasma density, and an anisotropic RD. The RD is seen to separate 
in time from the shockfront as expected. 


A. Formation of the ‘foreshock bubble’ 

The interaction of an RD with shock reflected ions was shown to lead to a compound 
structure referred to as the foreshock bubble (FB) which consists of a fast magnetosonic 
shock attached to a very hot, low density cavity behind it [25, 26]. The solution found here 
differs from the example of FB in Ref. [25] in at least two important aspects. First, the 
structure in Fig. 2 is more consistent with a localized front such as a shocklet than a shock 
since the plasma state behind it goes back to the solar wind conditions as evidenced from 
the ID cut of plasma variables (e.g., Figure 2d, 2j). Second, the cavity is seen to remain 
much smaller than the extent of the shockfront in y and does not grow much in time (Figures 
2a-2c). In contrast, in the case of FB in Ref. [25] the cavity covers almost the entire extent 
of the area behind the shock. Given these differences, we have labeled this structure as ‘FB’ 
to distingish it from the previously reported FBs. 

As it turns out, there is a spectrum of solutions between the example found here and 


that reported in Ref. [25, 26]. The resulting compound structure can take on several forms 
depending on the thickness of the RD, the IMF direction, among others. Particularly impor- 
tant is the level of density depression in the cavity formed near the shockfront. The larger 
the density depression, the hotter the temperature within the cavity. In moderate density 
depressions such as those in Figure 2 (> 0.1IV o ), the cavity does not increase in size signif- 
icantly. Note, however, that the size of the shockfront here in both x and y directions can 
still reach several hundred ion inertial lengths as seen in Figure 2. For more extreme density 
depressions approaching ~ 0.05-/V o , the shockfront and the cavity become one structure and 
can grow to very large scales even reaching the size of the magnetosphere if given sufficient 
time. At least this is as is observed in hybrid simulations (not shown). However, one has 
to be cautious regarding simulations of such structures in the hybrid code. Standard hybrid 
codes are not ideally suited for treatment of structures when the density depression reaches 
the density floor set in the simulations. And we have found a tendency to see weaker cavities 
with use of larger number of particles per cell and/or higher resolutions. We will examine 
formation of foreshock bubble using fully kinetic simulations elsewhere. An observational 
signature of such a structure, assuming it can be realized in nature, would be a large region 
of the solar wind populated with very low density and hot plasma. 

B. Characteristics of turbulence and its time evolution 

As is evident from Figure la, the quasi-parallel magnetosheath appears more turbu- 
lent than the quasi-perpendicular magnetosheath. And after the IMF change, the new 
quasi-parallel magnetosheath in Figures la-c appears more turbulent than the new quasi- 
perpendicular magnetosheath. This is confirmed by calculation of the magnetic field fluctu- 
ations which are about a factor of 4 larger in the Q\\ regions. This higher turbulence level is 
also confirmed by additional diagnostics in Figure 3. Figures 3a-b show the kurtosis (fourth 
standardized moment) of magnetic held increments 5B(x, s ) = B(x + s) — B(x) as a func- 
tion of the separation length s and the magnetic held spectrum for the foreshock, Q|| and 
Q± magnetosheath. Kurtosis, or hatness, is a measure of deviations from Gaussian. The 
hatness of Gaussian is exactly three. The Qy magnetosheath shows clear turbulent features, 
i.e. , an energy spectrum with close to —5/3 range and strong nonGaussianity (i.e., a very 
large kurtosis which increases as the length scale decreases). Interstingly enough, the kur- 
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tosis values are in good agreement with in situ measurements of Q\\ magnetosheath in cases 
where evidence of turbulent reconnection was reported [13]. In contrast, the kurtosis of the 
foreshock and the Q± magnetosheath are much smaller. In particular, the kurtosis of Q± is 
almost flat with values between 4 to 5, showing no or almost no signature of intermittency. 
Intermittency is a signature of the generation of coherent structures. As we will demonstrate 
shortly, turbulence generates coherent structures in form of reconnecting current sheets and 
magnetic islands in the Q\\ but not in the Q± magnetosheath, consistent with the kurtosis 
diagnostics. Figure 3c shows the comparison of the ion energy distribution function in the 
foreshock, Q||, and Q±. The foreshock shows the presence of counter streaming population 
as expected. The magnetosheath shows power-law distribution functions but with different 
power-law indices for the Q|| and the Q± cases. The thermal energy and the high energy 
particle production are significantly higher in the Q\\ magnetosheath. This is not surpris- 
ing since higher levels of turbulence in the Q\\ region can lead to additional heating and 
acceleration of particles. Another mechanism for acceleration to high energies is the Fermi 
process which requires the particles to go back and forth across the shock several times. 
Since particles can move along the magnetic field more readily than across it, one might 
expect more energetic particles in the Q\\ magnetosheath. The relative importance between 
these competing processes is beyond the scope of the present work and will be explored 
elsewhere. 

The power-law indices found here are not universal and can vary under different mag- 
netosheath conditions. Previous ID hybrid simulations of quasi-parallel shocks have also 
shown power- laws over a similar range [27, 28] but they show flattening of the particle spec- 
trum at higher energies where we do not have sufficient statistics. This flattening of the 
particle spectrum is also seen in 2D hybrid simulations of perpendicular shocks that include 
pre-existing, large-scale magnetic-field turbulence in the upstream region [29]. The details 
of the turbulence are clearly very different in our simulations as compared to ID simulations 
where reconnection and magnetic islands are excluded. Direct comparison of particle ener- 
gization in local shock simulations versus our global simulations will be presented elsewhere. 

In the present study, the initial solar wind ions are taken to have a Maxwellian distribution 
function. However, nonthermal ion distribution functions are commonly observed in the 
solar wind. The presence of nonthermals can affect the instabilities in the foreshock which 
in turn can affect the details of turbulence and associated particle energization. This is an 
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important extension of the present work that we are currently working on. 

Next, we compare the evolution of the turbulence level as measured by the root mean 
square of the magnetic held fluctuation amplitude B rms = ( < B >)"/n) which 

is normalized to B a . Here < B > is the average magnetic held in the region where B rms is 
calculated. We use run 2 where the IMF direction remains at 10° during the entire length of 
the simulation. We have run the simulation to longer times than previous hybrid simulations 
so that we can explore the asymptotic behavior of the turbulence. The result is shown in 
Figure 4a. Although the foreshock turbulence is a major source of turbulence downstream of 
the quasi-parallel shocks, the correlation between the turbulence level in these two regions 
is more complex as evident in Figure 4a. While on average the increase or decrease in 
the foreshock activity leads to a corresponding increase or decrease in the turbulence levels 
downstream, there is a temporal lag between changes in the foreshock and its transmission 
to magnetosheath. For example, there is no discernible increase in the foreshock turbulence 
levels from time of 200 to 500, whereas the Q\\ turbulence level has increased by a factor 
of nearly 2.4. However, in time the foreshock turbulence level starts to drop off, and the 
turbulence level downstream of the shock also starts to decrease. 

As another way to see this, and also to get a handle on the timescale for decay of 
turbulence levels once the geometry changes from quasi-parallel to quasi-perpendicular (Fig. 
4b), we examine time evolution of B rms for run 3. This run, which is a higher resolution 
(0.5di) version of run 1, has identical parameters as run 2 except that the IMF direction is 
changed from —45° to 10° after f l ci t — 90. We used this run as part of our convergence study 
to ensure the results are not affected by cell resolution or number of particles. It is seen 
that the change to quasi-perpendicular geometry results in a rapid drop in the turbulence 
fluctuation levels to about B rms ~ 0.3 which is similar to that in the Q± magnetosheath 
in Figure 4a. It is interesting to note that the mean rate of solar wind IMF directional 
discontinuities is about 0.5 — 2.5 hours [30] or about Q a t ~ 150 — 750. This is competitive 
with decay times of turbulence in the magnetosheath and provides a natural low frequency 
“stirring” mechanism for the magnetosheath turbulence. 

The significantly enhanced level of turbulence in the Q\\ magnetosheath has implications 
for the nature of turbulence as we now demonstrate. 
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C. Reconnection in the magnetosheath 


In 2D, the presence of magnetic islands formed through reconnection can be demonstrated 
by plotting the contours of the out-of-the-plane vector potential or tracing of magnetic field 
lines. An effective way to visualize field lines is through the line-integral-convolution (LIC) 
technique (see the Appendix) which we first introduced in plasma physics in our study of 
shear driven turbulence [3]. Figure 5 shows the magnetic field lines colored by V x for a 
zoomed in area of the quasi-parallel shock for run 1 at D ci — 694. The magnetic held 
intensity plot of a larger region of this simulation was shown in Figure lc. Figure 5 clearly 
shows a wide range of magnetic islands up to ~ 30 di in diameter. Some of the islands show 
evidence of outflowing jets (not shown). 

These magnetic islands are observed in the quasi-parallel but not in the quasi- 
perpendicular magnetosheath (not shown). The absence of magnetic islands in the Q± 
magnetosheath is the reason for differences in the kurtosis in the Q|| and Q± regions shown 
earlier (Fig. 3b). Also note that the islands can exist even in the vicinity of the Qy region 
of the bow shock. Due to the small size of these magnetic islands/flux ropes, they are not 
expected to be force free. 

These results are consistent with Cluster observations of magnetosheath reconnection 
events [12, 13] which also were only present in the quasi-parallel magnetosheath. However, 
they did not look for magnetic islands in the observations. Detection of small magnetic 
islands is challenging but given that some of the islands grow to as large as ~ 30 — 40 d i: it 
may be possible to find such islands in the data. 

We have found the reconnecting current sheets and magnetic islands to be a common 
feature of quasi-parallel shocks. Since the size of the magnetosheath is larger for larger D p , 
the number of islands and their interactions increases at larger D p . Figure 6 shows the 
density and the LIC of magnetic held colored by the strength of B for run 4 ( D p = 300). 
The magnetosheath is found to be almost completely volume filled with magnetic islands in 
this case. 


12 


D. Enhanced heating of Q y magnetosheath 


The enhanced turbulence level in Q|| magnetosheath also leads to stronger heating as 
compared to Q± magnetosheath as evident in Figure 7. We have found evidence for this 
in the observations. Figure 8 shows the statistical map of THEMIS measurements between 
October 2007 and October 2013 when the probes occupied the magnetosheath. The quan- 
tity binned here is the mean value of the total ion temperatures during a Parker-Spiral 
IMF ( B x > 0.2|F>| and B y < 0.2|F?| or B x < 0.2|R| and B y > 0.2|R|), which is presented 
in the Magnetosheath Interplanetary Medium (MIPM) reference frame. The MIPM frame 
uses upstream measurements to arrange within a normalised model magnetosheath thus ac- 
counting for boundary motion and planetary aberration. For a detailed description of the 
MIPM frame and subsequent data processing we refer the reader to access the methodology 
paper which also compares the statistical mapping results against global MHD simulations 
produced by the CCMC BATS-R-US code[31]. To briefly summarise, the MIPM frame ar- 
ranges points as a function of the upstream IMF. Therefore, statistically it can be closely 
compared to a aberrated GSE system during parker-spiral IMF orientation were Ymipm > 0 
represents a quasi-perpendicular shock and Ymipm < 0 quasi-parallel. The physical dimen- 
sion of each bin is 0.5 x 0.5 Earth radii and the bin number density is typically a hundred 
data points per bin. Each datapoint is calculated from the mean averaged 3-minute in- 
tervals of THEMIS data. Currently the data is binned using a 20 minute sliding averaged 
OMNI data around each 3 minute THEMIS data interval. The points during which IMF 
cone angle varies more than 20° degrees are removed. It is clear that the dayside MSH is 
hotter due to the shock compression and the hot plasma region extends more tail-ward at 
the quasi-parallel shock side of the magnetosheath, possibly due to kinetic processes. The 
quantitative calculation as a function of MIPM-frame zenith angle shows that the magne- 
tosheath downstream the quasi-parallel bow shock is on average about 10-15 percent hotter. 
A more detailed temperature study using THEMIS data for various solar wind conditions 
and quantitative comparisons between dawn and dusk sectors will be discussed elsewhere. 
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E. Global fully kinetic simulation 


Hybrid codes do not include electron kinetic effects. As such, magnetic reconnection in 
hybrid simulations occurs due to either numerical or imposed resistivity. To confirm that 
the magnetosheath reconnection events observed in our hybrid simulations are not due to 
numerical effects, we have also conducted a 2D global fully kinetic simulation for the radial 
IMF condition. The detailed comparison of this run with an equivalent hybrid simulation 
will be presented elsewhere. Here, our interest is only to establish that reconnection and 
magnetic islands also occur in the global fully kinetic simulation. 

The simulation was performed using a general-purpose plasma simulation code VPIC [32]. 
The simulation is of size 2000 x 1500 cl t with 56544 x 41984 cells. The IMF magnetic held is 
in the x direction and the dipole is oriented at an angle of 80° with respect to the IMF. The 
incoming solar wind plasma has parameters Ma — 8, T) = T e , (3 = 87r(T e + Ti)n 0 /Bl — 1, 
uj pe /(jj ce — 3. The solar wind is normal to the dipole. The ion-to-electron mass ratio is 
mi/m e — 50. The strength of the dipole corresponds to stand-off distance D p — 50 di. 
Reflecting boundary conditions for particles are prescribed on a circle of radius R = 0.3 D p 
around the geometrical center of the dipole. The dipole is created by prescribing out-of- 
plane current in two circular regions of the simulation domain or radius 0.1 R separated 
by distance 0.25A. The initial conditions correspond to a uniform plasma with solar wind 
parameters. In order to minimize the perturbations created by the initial conditions, the 
strength of the currents creating the dipole is ramped up in time as j — j 0 [l — exp(— t/r)], 
where r = SOD” 1 . 

The result is illustrated in Figure 9 where only a segment of the simulation around 
the quasi-parallel region is shown. The LIC of magnetic field is shown colored by ion 
density. As in the hybrid example, many magnetic islands are observed in the quasi-parallel 
magnetosheath. Our closer examination of the magnetic field structures (not shown) shows 
the hierarchy of islands extending down to the electron kinetic scales. However, the largest 
islands have similar size (~ 30 di) as in the hybrid simulations. 
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IV. TRANSMISSION OF THE FORESHOCK TURBULENCE 


It is well known that the relative streaming of ions reflected from the bow shock and the 
incoming solar wind leads to the generation of large amplitude ultralow frequency (ULF) 
waves [6]. As these waves get convected back towards the shock, they can steepen, giving 
rise to shocklets and short large-amplitude magnetic structures (SLAMS). Here we examine 
the evolution of such structures, their interaction with the bow shock, and their role in 
magnetosheath turbulence. For simplicity, here we refer to the steepened structures as 
wavefronts or shocklets. As seen from Figure 10, the wavefront-magnetosphere interaction 
falls into four stages. The IMF is nearly radial in this case (run 5). The part of the wavefront 
hitting the bow shock first gets compressed to shorter scales (Fig. 10 a,c,d), causes a rippling 
of the shockfront as well as local spike in the temperature (Fig. 10b). As the wavefront gets 
through the shock (Fig. lOe), its amplitude is amplified (Fig. lOg). The waves persist as 
wavefronts in the magnetosheath for tens of ion cyclotron times (Fig. 10i-l) but then they 
start to break up and dissipate their energy. In this way, turbulence upstream gets amplified 
and contributes to the turbulence in the magnetosheath. 

Interestingly enough, large amplitude waves have also been observed downstream of the 
quasi-parallel shocks at Venus. A recent study compared detailed wave diagnostics in both 
upstream and downstream regions of the Venusian bow shock and concluded that the down- 
stream waves are transmitted ULF waves [33]. 

In cases where the steepened fronts are nearly parallel to the quasi-parallel region of the 
bow shock surface, they can significantly modify the shock structure and the magnetosheath 
thickness. This is shown in Fig. 11, from run 6, which shows the time evolution of density and 
B y . The propagation of one wavefront, labeled as WF1, is marked in time. The steepened 
front is evident in Figures 11a and lie. Note that the steepened front has many properties 
similar to a shock, exhibiting localized density and magnetic pulses/enhancements over the 
background values. By the time — 47.2, the front has gained strength (Figures lib and 
Ilf) and creates a large density jump ahead of the shock (Figure lib). The thickness of the 
magnetosheath is much thinner in this region as compared to the upper region away from 
the shocklets. A rather large magnetic island of about 10 di in width is evident as a circular 
loop in the density panels. WF1 grows in scale (Figures 11c and llg) and eventually gets 
convected into the shock and moved downstream (Figures lid and llh). This process can 
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repeat for several cycles. This is evident in Figures lid and llh which show the formation 
of a second wavefront (WF2) which will go through the same process as WF1. Figure 12 
shows a blow up of B y and an associated ID cut across the wavefront at f l C jt — 47.2. The 
ID cut shows the amplification of the wavefront amplitude as it crosses the shock into the 
magnetosheath. 

Figure 13a-b illustrates the dynamics between the shocklet and the bow shock. The LIC 
of magnetic field colored by density and LIC of ion velocity colored by V x are shown for a 
zoomed in region of Figure 11c. The shocklet presents an obstacle to the flow, generating 
local spikes in the density (Fig. 13a) and the magnetic field, and decelerating the flow 
(Fig. 13b). The region between the shocklet and the magnetosheath can be quite turbulent. 
Several vortices and magnetic islands are clearly evident in this in-between region, similar 
to the magnetosheath. Note also the presence of a relatively large magnetic island spanning 
about 10 di along the magnetopause surface. The vortices are apparent here since the flow 
speed is nearly zero in the shown regions. The crossing of such shocklet /bow shock structures 
may be mistaken as multiple shock crossings in the data, especially since there are regions 
with conditions similar to the solar wind in between them. 

Another effect of ion foreshock that is perhaps not widely appreciated is that it reduces 
the magnetosonic Mach number as evident from the ID cut in Figure 14. As such, the ion 
foreshock helps in the shock dissipation process and affecting the slowing of the plasma. The 
presence of the jets can also been as streaks of high magnetosonic Mach number within the 
magnetosheath (Figure 14). The modification of M s due to the ion foreshock is consistent 
with an observational study where significant evolution of interplanetary shock parameters 
were observed during its motion along the bow shock [34]. 


V. ANOMALOUS FLOWS 

In this section we demonstrate the effects of turbulence on the magnetosheath flow. There 
are a number of diagnostics that are useful for visualizing flow. These include streamlines, 
pathlines, and streaklines. Streamlines represent instantaneous patterns of the flow. A true 
vortical pattern on scales larger than kinetic scales would wrap the magnetic field. However, 
the presence of such a vortex can be masked, using streamlines, in the presence of a large 
background flow as is the case in the magnetosheath. Instead, we use streaklines to study the 
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mechanisms for formation of sunward flow and vortices in the magnetosheath. Streaklines 
provide a time history of the flow and can reveal the presence of vortical pattern even in 
case of large background flow. Each streakline is created by tagging and following the trace 
of fluid elements passing through a particular spatial point. Streaklines can be thought of 
the pattern traced by dye if one were to steadily inject dye into the fluid at a fixed point. 

A. High speed jets 

An interesting feature evident in Figure 5 is the presence of two ion jets, shown as red, 
that penetrate deep inside the magnetosheath. The magnetic held lines are seen to remain 
nearly radial as in the solar wind along these jets, embedded in an otherwise sea of magnetic 
islands. In this section, we examine these jets in more detail. Figure 15 shows plots of 
the dynamical pressure over a larger region of the simulation than in Figure 5. Presence 
of many jets is evident, seen as streaks of high dynamical pressure, with some reaching 
the magnetopause. As the jets slow down in the magnetosheath, they form a bow wave 
(Fig. 15) which can be as large as several hundred d*. The interaction of the jets/bow wave 
with the magnetopause causes strong perturbations, giving rise to surface waves, and/or 
inner-magnetospheric compressional waves [11, 35, 36]. And in this case, this interaction 
has triggered flux transfer events (Fig. 15b), although the IMF is 10°. The jets also lead 
to significant inhomogeneities in the magnetosheath and serve as an additional driver of 
turbulence. Figure 15b shows the LIC of B for an area zoomed in around the large bow 
wave in Fig. 15a. This leads to another source of reconnection and generation of magnetic 
islands in the magnetosheath. 

In time, the jets get convected towards the tail while new jets are formed. Figure 16 shows 
the properties of the jets, which strongly deviate from the surrounding magnetosheath. In 
particular, jets exhibit enhanced magnetic held strength (Fig. 16b), lower temperature 
(Fig. 16c), and higher density (Fig. 16d) than the surrounding magnetosheath. Note also 
the strong kinking of the magnetic held (Fig. 16b) along the jets. In each case, we have found 
the wavelength of the kinked jet to be in good agreement with the wavelength of Kelvin- 
Helmholtz instability based on local measurements of velocity shear and the thickness of the 
shear layer. 

Observations have revealed the presence of transient enhancements of plasma how within 
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the magnetosheath characterized by various measures such as high dynamic pressure and 
kinetic energy density [9, 10, 35, 37-40]. Hietala et al. [9] reported high velocity jets during 
nearly radial conditions (IMF angle less than 20°) for Ma — 12. They found the plasma 
density to be compressed along the jets with density depletions around them, the velocity 
of the jets to be close to the upstream value, and the width of them on the order of 50 to 
100 ion di . These features are fully consistent with the observed jets here. The velocities of 
the jets shown in Figure 5 are close to the solar wind speed of 8 Va, their widths are in the 
range of 20 — 60 d,, and there is density compression on the order of a factor of 6 along the 
jets (Fig. 16d) and large density depletions around them, with densities dropping below the 
solar wind density in some cases. The size and the level of the depletion area varies in time. 
The ion temperature (Fig. 16c) shows that the temperature does not change significantly 
along the jets. This is also consistent with properties of the observed jets [9, 10]. This may 
also explain previous observations [41] which had found cases where the solar wind passes 
through the shock layer without significant heating. 

The origin of these jets has been controversial. Hietala et al. [9] attributed their formation 
to ripples along the shock surface. Since the shock mainly decelerates the component of the 
upstream velocity parallel to the shock normal, ripples along the shock cause local changes 
to the angle between these two vectors. For large angles, the shock speed would remain 
close to its upstream value. While others have suggested the origin of the jets to be due to 
interaction of the magnetosphere with solar wind discontinuities (e.g., Ref. [42]) and others 
(e.g., [39]) questioned the ripple-based mechanism. In contrast, a more recent study [10] 
using 4 years of subsolar magnetosheath observations concluded that 97% of the observed 
jets are consistent with the ripple idea. This finding is consistent with our simulation results 
where we find the bow shock surface to be highly rippled due to foreshock turbulence, and 
we find the jets even in cases where we keep the IMF direction fixed in time and there are 
no solar wind discontinuities. 

B. Sunward Flows 

In the absence of turbulence, the downstream flow must be directed away from the bow 
shock at all points as seen in global MHD simulations. However, THEMIS observations 
during radial IMF conditions have shown evidence of sunward flows. Such flows are also 
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seen in our simulations. Figure 17 shows V lx for the same time in the simulation as in Figure 
16. The purple areas are regions in the magnetosheath with negative V lx with peaks around 
2Va- Note that there are several regions along the bow shock where such sunward flows are 
present. 

Several different mechanisms have been proposed to explain this anomalous flow. Shue 
et al. [11] attributed the origin of the sunward flow to the interaction of a fast anti-sunward 
flow of unknown origin with the magnetopause, which in turn causes the magnetopause to 
rebound, creating a sunward motion. As shown in Figure 16, the generation of fast anti- 
sunward flow is due to the rippling of the bow shock. Although some of the jets interact with 
the magnetopause, many terminate at significant differences away from the magnetopause 
(e.g., Figure 16a) and yet there are sunward flows in the vicinity of such jets (Fig. 17). This 
indicates that alternative mechanisms must be in play for the generation of sunward flow 
than that proposed by Shue et al. [11], The simulations enable us to study the generation 
mechanism of sunward flow in much more detail than is possible in spacecraft measurements. 

We have seen that the jets snake around the highly inhomogeneous magnetosheath, some 
reaching the magnetopause while others are terminated earlier. As we will see, some of the 
jets encounter local magnetic held that is nearly perpendicular to them. This causes them 
to get deflected since it is difficult for plasma to traverse perpendicular to the magnetic 
held. This complex dynamics has a direct bearing on the formation of sunward how. We 
use streaklines to follow the evolution of the how. First, we examined the role of jets that 
reach the magnetopause in causing the sunward how. Figure 18a shows the streaklines 
shortly after they were seeded in two regions in the foreshock. These two groups of seeds 
were carefully selected so as to follow both the plasma motion inside the two large jets that 
reach the magnetopause (painted in white) as well as the plasma motion aroung the jets 
and their bow waves (painted in blue). Several points are immediately clear in Figure 18. 
First, the sunward how is not caused by reflection of the jets. Second, there is no evidence 
that magnetopause motion is creating a large scale sunward how. Third, in Figure 18d, one 
of the largest jets has started to wrap although it has not formed a complete vortex. Any 
wrapping clearly created a sunward how locally. Fourth, it is evident in Figures 18b-d that 
the large bow wave pushes on the plasma around it, seen as blue streaks, as it propagates 
downward. This is another way that sunward hows can be generated. 

Next, we seeded the held lines to examine the origin of the sunward hows between the jets 
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that extend into the solar wind as were seen in Figure 17. The evolution of the streaklines 
is shown in Figure 19. Here again we marked the seeds in two colors, those that interact 
with the two jets marked as J\ and J 2 , marked in white and those that are from the areas 
in the vicinity but outside of these jets, marked in blue. The streaklines at the mouth of the 
opening in the region between ,J\ and J 2 are seen to bend sunward as the blue streaklines 
penetrate this area (Fig. 19c-f). The blue streaklines represent plasma that was pushed 
downward due to the pressure of the jets on the ambient plasma. In this way, the jet energy 
causes deflection of the ambient sunward, rather than the plasma within the jets being 
deflected. 

As another check on the formation mechanism of sunward flow, we show in Figure 20 a 
zoomed-in view of the region containing jets marked as J\ and J 2 in Figures 19d and 19e. 
The top two panels show the LIC of magnetic held colored by ion temperature. The jets 
are seen to have a temperature similar to the solar wind whereas the area between the jets, 
marked as R2, is much hotter. The magnetic held lines clearly show that the surface of 
the bow shock is highly corrugated. A few sample magnetic held lines are plotted in white, 
while density contours are plotted in black. Note the presence of a band of nearly laminar 
magnetic held that is nearly perpendicular to the jets as they enter the magnetosheath, 
marked as Rl. As the nearly parallel propagating jets encounter this band of magnetic held, 
they cannot penetrate further into the magnetosheath as it is difficult for plasma to move 
across the magnetic held. The upper jet is seen to snake around in the magnetosheath and 
gets dehected when reaching this band of magnetic held. The lower jet also winds around 
but does not cross this band of magnetic held. The plasma in the region between the two 
jets -Rl gets squeezed and follows the magnetic held lines downward. This plasma pushes 
outward the shocked solar wind between the two jets in region R2, causing sunward how as 
evident from negative Vi x hows seen in Figure 20 c-d. This picture is consistent with the 
development of the how seen in Figure 19. The blue streaklines in Figure 19 represent the 
how that follows the band of magnetic held lines Rl and enters the region between the two 
jets R2, creating sunward how. 

These considerations have led us to the following picture for formation of sunward hows. 
The jets, which have high dynamical pressure, “stir” the magnetosheath, creating large scale 
disturbances and pushing the ambient lower density plasma and the embedded magnetic 
held. The plasma in the sunward how regions is not the plasma in the jets that was dehected, 
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rather it is the plasma that was already there that is accelerated by thermal and dynamical 
pressure that develops because of the jets. In other words, the energy of the jets is deflected 
into sunward flow rather than the deflection of the jets themselves. This energy deflection 
can give rise to sunward flow that extends tens and even hundreds of ion inertial length into 
the solar wind, well past the nominal position of the bow shock (Fig. 17). 

C. Large scale vortices 

The presence of jets and sunward flows in the Q|| magnetosheath can drive velocity shears 
which can lead to formation of large vortices. Figure 21 (run 7) shows the formation of a 
large scale vortex which occurs in two stages. A structure akin to the Rayleigh- Taylor finger 
is formed (Fig. 21a) which then gets caught into the velocity shear in the magnetosheath 
(Fig. 21b), forming a large scale vortex (20c). The magnetic held lines are wrapped around, 
similar to a vortex generated by Kelvin-Helmholtz (KH) instability. Figure 22 illustrates 
the generation mechanism of the vortex using streaklines. A band of streaklines above the 
nominal stagnation line how is seen to get deflected sharply and turn downward (pink) (Fig. 
22a). This is due to the dehection of the energy of a jet which pushes the ambient plasma 
downward (not shown) similar to that for run 1 in Figure 19. However, unlike in Figure 
19, here this plasma gets caught in the velocity shear in the magnetosheath (Fig. 22b) and 
turns into a vortex (Fig. 22c). The color coding of streakline is also useful for tracking the 
movement of the plasma in different regions of the magnetosheath as seen in Fig. 22c. The 
structure in Figure 21a has the appearance of a Rayleigh-Taylor finger. The plasma in the 
jet is denser than the ambient plasma which in the presence of an accelerating force that 
would mimic gravity could give rise to the Rayleigh-Taylor instability. However, here the 
formation of this structure appears to be due to dynamics, i.e. , deflection of the jet energy 
as it encounters a nearly perpendicular magnetic held. Although KH driven vortices are 
known to form on the hanks of the magnetopause (e.g., Ref. [18] and references therein) due 
to the velocity shear between the solar wind and the magnetopause, we believe this is the 
hrst report of a vortex driven due to turbulence in the magnetosheath. 
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VI. DISCUSSION AND SUMMARY 


The fundamental processes of magnetic reconnection, shocks and turbulence have been 
demonstrated to be more intimately connected than previously realized. Laminar reconnect- 
ing current sheets can develop turbulence [43] ; studies of turbulence show that formation of 
reconnection current sheets may be a generic feature of turbulence in magnetized plasmas 
(e.g., [1-5] and references therein). The link between shocks and the other two processes 
has, until this study, been less clear. For example, one early model of reconnection, the 
Petschek model, required four standing slow shocks as part of the reconnection exhaust [44]. 
However, the applicability of this model in collisionless plasmas remains in question [45]. 
Turbulence does not lead to shock formation, except possibly in shear driven turbulence [3], 
but even there turbulence is not the cause of shock formation. 

Shocks of all types in collisionless plasmas share the property of generating reflected dis- 
turbances of electrons and especially ions that can get back upstream into the unshocked 
medium [6, 46]. The angle of upstream magnetic held to the shock normal guides the re- 
fected particles and is one of several parameters controlling the spatial volume to which they 
may reach before being swept back through the shock layer. In addition the shock layer sub- 
structure may provide for the instability of the transmitted distributions, as is particularly 
well known for quasi-perpendicular shocks. The turbulence on the downstream side of the 
shocks reflects the geometry of B to the normal, with the quasi-parallel shocks possessing 
downstream fluctuations comparable to the mean field and the quasi-perpendicular down- 
stream regions quieter. We have demonstrated for the first time that turbulence associated 
with the Q y shocks lead to generation of magnetic islands, with some islands even forming at 
the shock surface. In some cases (e.g., Fig. 13) the islands can even form ahead of the shock. 
We have also demonstrated that reflected upstream disturbances when converted back and 
transmitted through the shock layers eventually dissipate there. Thus both ’’waves” and 
reconnection play a role in providing dissipation downstream of the quasi-parallel shocks. 
The question of the relative importance of these two processes for dissipation is beyond the 
scope of this paper and will be addressed elsewhere. 

We have thus forged a link between shocks, turbulence and magnetic reconnection in 
collisionless plasmas. In so doing we have provided a theoretical argument that explains 
the location of reported reconnection events in the magnetosheath [12, 13]. Previous studies 
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have not reported the incidence of reconnection in simulations of fast shocks or global hybrid 
simulations possibly because (i) many fast shock studies were in 1-D where reconnection 
cannot occur; (ii) the very limited number of 2-D quasi-parallel shock simulations apparently 
did not check for incidence of reconnection; and (iii) the low number of particles per cell and 
or low grid resolution of these early simulations (e.g., [25]) may have masked the detection of 
magnetic islands. Our findings here show that with increasing turbulence levels the system 
exhibits greater structural organization in form of coherent structures such as flux ropes, 
vortices and fast jets. Whether such reorganization into coherent structures is a consequence 
of a general governing principle of all systems far from equilibrium that are driven remains 
an intriguing speculation. 

Another interesting result here is in regards to the large scale consequences of the bow 
shock turbulence in the magnetosphere. This turbulence can, under certain conditions, 
enable formation of large scale jets that can interact with the magnetopause and lead to 
surface waves on it, or inner magnetospheric compressional waves [11, 35, 36]. The jets can 
also lead to formation of sunward flows and large scale vortices downstream of the shock. 

Given that we found magnetic islands to be a common feature of quasi-parallel shocks 
with B rms ~ 1, search for such structures in the data would be of interest. Ordinarily flux 
transfer events are conceptually associated with southward IMF ; our simulations have shown 
examples where flux transfer events can occur with northward IMF along the magnetopause 
(Fig. 15b) as well as inside the magnetosheath (e.g., Fig. 6b). While data identified flux 
transfer events can reach scales of hundreds of di , our results demonstrate that flux ropes 
due to shock ‘turbulence’ appear to be smaller and at most reach tens of di in scale. 

We have already started the extension of the present work to 3D. Our preliminary re- 
sults indicate that the basic physical effects uncovered here, such as reconnection in the 
magnetosheath, remain intact in 3D. However, 3D introduces additional complexity such as 
fanning of the magnetic field lines in the ion foreshock and formation of Kelvin-Helmholtz 
instability along the flanks. These results will be reported elsewhere. 

APPENDIX A: LINE INTEGRAL CONVOLUTION TECHNIQUE 

The standard way to visualize streamlines of a vector field is to seed some points and 
integrate to trace curves that are instantaneously tangent to the velocity vector. Although 
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useful, this approach can be cumbersome when interactively exploring a dataset. For ex- 
ample, the technique is inherently local and unless feature locations are know apriori it 
is difficult to select an effective set of seed points. An alternative technique called “line 
integral convolution or LIC, which convolves noise with a vector field producing streaking 


streamlines over the entire computational domain is found in one step without the need to 
explicitly specify a set of seed points. 

The technique was initially developed for use on images [47], but has since been extended 
to arbitrary surfaces[48], and to volumes [49]. The algorithm has been ported to the GPU[50] 


nique and have made it available through the latest release of ParaView[52], an open source, 
cross platform, tool for parallel interactive visualization of large datasets. In addition to tra- 
ditional visualization algorithms ParaView provides an MPI based data-parallel GP-GPU 


rithm designed for interactive data exploration. ParaView’s surface LIC can be used on 
massively parallel supercomputers without GPUs. It is described in detail in ParaView’s 
online documentation [53]. To aid the reader in the interpretation of the figures presented 
here a brief review of the LIC algorithm is presented here. 

The LIC algorithm of a vector field defined on an image, I(x,y ) is given by the following 
integral over streamline arcs computed from the center of each pixel location, x, y. 


where L is the integration length, N generates the noise value at a given location, Si is 
a position on a streamline arc centered at x, y, k(i ) is an appropriate convolution kernel, 
and I(x,y ) is the image pixel at x,y. In practice, streamline arcs are computed using an 
RK method over a fixed number of steps, and L can be a constant for all pixels in I, or it 
may be a function of the local vector field. When L is constant the resulting visualization 
has a uniform look. When L varies as a function of the vector local field the visualization 
accurately shows relative strengths of flow features. Because the LIC produces a dense 
representation of the flow field, features of interest can be quickly identified. LIC is often 
used in conjunction with scalar pseudocoloring. In this case a specialized shader is used to 


patterns that follow vector field tangents, has the advantage that a very detailed view of the 


and data-parallel implementations have been developed [51]. We have extended this tech- 


surface LIC algorithm which includes a number of enhancements to the basic LIC algo- 
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combine the colors with the gray scale LIC image. Optimized implementations, especially 
data-parallel or GPU based, are very fast making the technique useful for interactive data 
exploration. 

An example of image LIC of electron stream lines colored by the magnitude of electron 
velocity is shown in figure 23a. The data is from a 2D fully kinetic simulation of asymmetric 
current sheet. The formation of a vortex associated with a magnetic island (not shown) is 
clearly evident. The surface LIC technique we use is similar to image based LIC except first 
vectors are projected onto the surface then into image space where an image LIC algorithm is 
used to compute LIC [48]. Lit, pseudocolored surface geometry, is combined with the image 
LIC to produce a realistic rendering of the surface geometry. To show the utility of surface 
LIC in analysis of 3D data, we have applied it to data from a 3D fully kinetic simulation 
of reconnection [43]. Figure 23b, using isocontour of density, shows the presence of multiple 
flux ropes as evident from their helical magnetic field lines. 


ACKNOWLEDGMENTS 

This work was partially supported by NASA (NNX12AD30G, Heliophysics Theory), NSF 
(0904734, 1104815), and DOE (DE-SC0004662). M. W. acknowledges support from NSF 
Grants No. AGS-1063439 and No. AGS-1156094 (SHINE), and from NASA Grant No. 
NNX11AJ44G. Simulations were performed on Pleiades provided by NASA’s HEC Program, 
and Blue Waters sustained-petascale computing project, which is supported by the NSF 
(OCI 07-25070) and the state of Illinois. We gratefully acknowledge useful discussions with 
J. Giacalone, J. Borovsky, and A. Retino. We also thank the referee for making useful 
suggestions that improved the paper. 


[1] W. H. Matthaeus and S. L. Lamkin, Phys. Fluids 29, 303 (1985). 

[2] S. Servidio, P. Dmitruk, A. Greco, P. Wan, S. Donato, P. A. Cassak, M. A. Shay, V. Carbone, 
and M. Matthaeus, Nonlinear Processes Geophys. 18, 675 (2011). 

[3] H. Karimabadi, V. Roytershteyn, M. Wan, W. H. Matthaeus, W. Daughton, P. Wu, 
M. Shay, B. Loring, J. Borovsky, E. Leonardis, et al., Physics of Plasmas 20, DOI: 

10.1063/1.4773205, 012303 (2013). 


25 


[4] H. Karimabadi, V. Roytershteyn, W. Daughton, and Y.-H. Lin, Space Science Review DOI 
10.1007 / si 1214-013-0021-7, 178 (2013). 

[5] H. Karimabadi and A. Lazarian, Physics of Plasmas 20, 112102 (2013). 

[6] D. Burgess and et al., Space Science Review 118, 205 (2005). 

[7] B. Lembege and J. Dawson, Phys. Fluids 30, 1767 (1987). 

[8] O. Alexandrova, Nonlin. Processes Geophys. 15, 95 (2008). 

[9] H. Hietala, T. V. Laitinen, K. Andreov, R. Vainio, A. Vaivads, M. Palmroth, T. I. Pulkkinen, 
H. E. J. Koskinen, E. A. Lucek, and H. Rme, Physical Review Letters 103, 245,001 (2009). 

[10] H. Hietala and F. Plaschke, Journal of Geophysical Research 118, 7237 (2013). 

[11] J.-H. Shue, J.-K. Chao, P. Song, J. P. McFadden, A. Suvorova, V. Angelopoulos, and 
F. Plaschke, Geophys. Res. Lett. 36, L18112 (2009). 

[12] A. Retino, D. Sundkvist, A. Vaivads, F. S. Mozer, M. Andr, and C. Owen, Nature Physics 3, 
235 (2007). 

[13] D. Sundkvist, A. Retino, A. Vaivads, and S. D. Bale, Phys. Rev. Lett. 99, 025004 (2007). 

[14] B. Sonnerup, H. Hasegawa, W.-L. Teh, and L.-N. Hau, Journal of Geophysical Research 111, 
A09204 (2006). 

[15] N. R. Hasegawa, H., M. Fujimoto, V. Sergeev, E. A. Lucek, H. Rme, and Y. Khotyaintsev, 
Journal of Geophysical Research 112, A11206 (2007). 

[16] C. Russell and R. C. Elphic, Geophysical Research Letters 6, 33 (1979). 

[17] H. Hasegawa, M. Fujimoto, T.-D. Phan, H. Reme, A. Balogh, M. W. Dunlop, C. Hashimoto, 
and R. TanDokoro, Nature 430, 755 (2004). 

[18] T. Nakamura, W. Daughton, H. Karimabadi, and S. Eriksson, Journal of Geophysical Research 
118, 1 (2013). 

[19] H. Karimabadi, H. X. Vu, D. Krauss-Varban, and Y. Omelchenko, Numerical Modeling of 
Space Plasma Flows: ASP Conference Series, Edited by N.V. Pogorelov and G.P. Zank 359, 
257 (2006). 

[20] H. Karimabadi, B. Loring, H. X. Vu, Y. Omelchenko, M. Tatineni, A. Majumdar, U. Ayachit, 
and B. Geveci, 5th International Conference of Numerical Modeling of Space Plasma Flows, 
Proceedings of a 5th international conference held at San Diego, California, USA 13-18 June 
2010. Edited by Nikolai V. Pogorelov. San Francisco: Astronomical Society of the Pacific p. 
281 (2011). 


26 


[21] H. Karimabadi, P. Oleary, T. Mahidhar, B. Loring, A. Majumdar, and B. Geveci, XSEDE13 
Extreme Science and Engineering Discovery Environment: Gateway to Discovery San Diego, 

CA, USA July 22 - 25, 2013 ACM 978-1-4503-2170-9/13/07, New York, NY, USA, 

257 (2013). 

[22] J. D. Scudder and S. Olbert, Journal of Geophysical Research 84, 6603 (1979). 

[23] E. Sittler and J. D. Scudder, Journal of Geophysical Research 85, 5131 (1980). 

[24] N. Omidi, X. Blanco-Cano, C. Russell, and H. Karimabadi, Adv. Space Res. 33, 11 (2004). 

[25] N. Omidi, J. Eastwood, and D. G. Sibeck, Journal of Geophysical Research 115, A06204 

( 2010 ). 

[26] D. Turner, N. Omidi, D. G. Sibeck, and V. Angelopoulos, Journal of Geophysical Research 
118, 1552 (2013). 

[27] J. Giacalone, D. Burgess, S. Schwartz, and D. Ellison, Geophys. Res. Lett. 19, 433 (1992). 

[28] J. Giacalone, Astrophysical Journal 609, 452 (2004). 

[29] J. Giacalone, Astrophysical Journal 628, L37 (2005). 

[30] L. Burlaga, Solar Physics 7, 54 (1969). 

[31] A. P. Dimmock and K. Nykyri, J. Geophys. Res. 118, 4963 (2013). 

[32] K. J. Bowers, B. J. Albright, L. Yin, B. Bergen, and T. J. T. Kwan, Physics of Plasmas 15, 
7 (2008). 

[33] L. Shan, L. Q., M. Wu, X. Gao, C. Huang, T. Zhang, and W. S., Journal of Geophysical 
Research DOI: 10.1002/2013JA019396 (2014). 

[34] L. Prech, Z. Nemecek, and J. Safrankova, Earth Planets Space 61, 607 (2009). 

[35] E. Amata, S. P. Savin, D. Ambrosino, Y. V. Bogdanova, M. F. Marcucci, S. Romanov, and 
A. Skalsky, Planet. Space Sci. 59, 482 (2011). 

[36] F. Plaschke, K.-H. Glassmeier, D. Sibeck, H. U. Auster, O. D. Constantinescu, V. Angelopou- 
los, and W. Magnes, Ann. Geophys. 27, 4521 (2009). 

[37] Z. Nemecek, K. S. M. T. S. D. J. Safrankova J., Prech L., L. Prech, M. T. Kokubun, S., and 
D. Sibeck, Geophys. Res. Lett. 25, 1273 (1998). 

[38] S. Savin, E. Amata, L. Zelenyi, V. Budaev, G. Consolini, R. Treumann, E. Lucek, 
J. Safrankova, Z. Nemecek, Y. Khotyaintsev, et ah, JETP Lett. 593, 87 (2008). 

[39] M. O. Archer and H. T. S., Ann. Geophys. 31, doi:10.5194/angeo-31-319-2013, 319 (2013). 

[40] F. Plaschke, H. Hietala, and V. Angelopoulos, Ann. Geophys. 31, 1877 (2013). 


27 


[41] J. T. Gosling, M. F. Thomsen, S. J. Bame, and C. T. Russell, J. Geophys. Res. 94, 10027 
(1989). 

[42] S. Savin, E. Amata, L. Zelenyi, V. Lutsenko, J. Safrankova, Z. Nemecek, N. Borodkova, 
J. Buechner, P. W. Daly, E. A. Kronberg, et al., Ann. Geophys. 30, 1 (2012). 

[43] W. Daughton, Y. Roytershteyn, H. Karimabadi, L. Yin, B. J. Albright, B. Bergen, and K. J. 
Bowers, Nature Physics 7, 539 (2011). 

[44] H. Petschek, Proceedings of the AAS-NASA Symposium held 2830 October 1963 at the God- 
dard Space Flight Center, Greenbelt, MD p. 425 (1964). 

[45] A. Le, J. Egedal, J. Ng, J. H. Karimabadi, Scudder, V. Roytershteyn, W. Daughton, and 
Y.-H. L. Liu, Physics of Plasmas 21, 012103 (2014). 

[46] H. Karimabadi, D. Krauss-Varban, and N. Omidi, Journal of Geophysical Research 100, 11957 
(1995). 

[47] B. Cabral and L. Leedom, SIGGRAPH ’93. doi:10. 1145/166117.166151, 263 (1993). 

[48] R. Laramee, B. Jobard, and H. Hauser, Visualization, 2003. VIS 2003. IEEE 
doi:10.1109/VISUAL. 2003. 1250364, 131 (2003). 

[49] V. Interrante and C. Grosch, NASA (Hampton, VA US) (1997). 

[50] D. Weiskopf, Springer- Verlag New York, Inc., Secaucus, NJ, USA (2006). 

[51] S. Bachthaler, M. Strengert, D. Weiskopf, and T. Ertl, In proceeding of: Eurograph- 
ics Symposium on Parallel Graphics and Visualization, EGPGV 2006, Braga, Portugal 
DOL10.2312/EGPGV /EGPGV06/075-082 (2006). 

[52] A. Henderson, Publisher: Kitware, Inc. (2005). 

[53] Paraview, http://paraview.org/Wiki/ParaView/Line_Integral_Convolution (2014). 


28 




tfia=694 


3500 


Y(d i ) 


1000 



300 


Xldp 


1800 300 


X ( d .) 1300 300 x(d . } 1800 


B 


7 










= 0 


FIG. 1: (a)-(c) Time evolution of the magnetosphere in run 1 where the IMF direction changes in 
time. The arrows indicate the direction of the magnetic held. A large foreshock bubble is evident 
in panel a). Details of this structure is shown in the next figure. The positions of the Qy and 
Q± regions move due to the change in the IMF direction. Note the higher level of fluctuations 
associated with the Q y magnetosheath. In panel c), streaks of enhanced magnetic held are observed 
in the magnetosheath. These are associated with the formation of high velocity jet that penetrate 
deep into the magnetosheath. 
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FIG. 2: Time evolution of the foreshock bubble in run 1 and ID cuts of density, total ion temper- 
ature, the x and y compoments of ion velocity 1 4 r, Viy > and the y-component and total magnetic 
field By , and B. The change in the IMF direction launches a rotational discontinuity (RD) which 
interacts with the ions reflected from the shock, giving rise to a shockfront and a density cavity 
in the solar wind. The red arrows indicate the direction of the magnetic field. Reflected ions in 
the foreshock penetrate through the shockfront, streaming along the new incoming IMF direction. 
Observational signature is sunward moving particles on field lines that do not point at the bow 
shock. The RD separates from the shockfront in time. 
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FIG. 3: Comparison of turbulence characteristics in the foreshock, Q y magnetosheath and the Q± 

magnetosheath at f l c i — 640 for run 1 shown in Figure 1. a) Omnidirectional energy (per unit 

mass) spectra of magnetic held. Only Qy magnetosheath shows a 5/3 power- law in this case, b) 

Kurtosis (fourth standardized moment) of magnetic held increments 5B(x, s ) = B(x + s) — B{x) 

as a function of the separation length s is shown. Kurtosis is a measure of nonGaussianity of the 

turbulence. An increase in kurtosis at smaller scales is an indicator of the presence of coherent 

structures in turbulence. Only the Q|| magnetosheath shows this behavior in this run. c) Ion 

distribution function as a function of energy normalized to upstream how energy. Foreshock shows 

the presence of two populations (counter streaming). Both the Qy and Q± magnetosheath show 

power-laws. We have used the following range of x and y for each region in panels a)-b): for 

foreshock x — [301 — 556 di], and y — [1301 — 1812 di], for the Qy, x — [681 — 808dj], and y — 

[1651 — 1906di], and for the Q± x — [1401 — 1656], and y — [2801 — 3056di]. In panel c), we 

have binned the energy, measured in the simulation frame, into 500 logarithmically based segments 

and chosen the following range of x and y for each region: for foreshock x — [400 — 500di], and 
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y — [1500 — 2000di], for the Qy, x — [700 — 800di], and y — [1500 — 2000 di], and for the Q± 


x m [1500 - 1800dj], and y = [2700 - 2900^]. 


a) 




FIG. 4: Time evolution of turbulence fluctuation levels as measured by B rms which is normalized 
to B 0 . a) As the magnetosphere forms (run 2), ions start to get reflected from the Q\\ segment of 
the bow shock. This in turn drives the upstream turbulence as evident in an initial rise in B rms 
until it reaches a peak of nearly 2. After the peak, B rms starts to weaken as the foreshock B rms 
also weakens. B rms in the Q± is driven by the temperature anisotropy and is more steady in 
time. The higher B rms in Q y as compared to the foreshock is not surprising since the turbulence 
convected from upstream through the shock gets amplified in the process, b) In run 3, the IMF 
direction changes in time, resulting in the switch of the geometry of an originally Q || to Q±. This 
switch is shown to lead to rather rapid drop in B rms and it settles down to B rms levels similar to 
those in the Q± in panel a). 
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FIG. 5: Formation of a hierarchy of magnetic islands due to Q\\ magnetosheath turbulence in a 
global hybrid (runl) simulation. This is a zoomed view of the Q || region of Figure lc. Here the LIC 
technique is used to visualize the magnetic field lines. The LIC is colored by Vj x . In addition to 
the magnetic islands, the presence of at least two jets inside the magnetosheath having velocities 
comparable to the solar wind speed are evident (red streaks in the magnetosheath). 
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FIG. 6: Formation of turbulence and associated magnetic islands for a run with D p — 300 (run 
4). We have examined the properties of turbulence for a range of parameters, including different 
dipole field strength. The formation of the magnetic islands seems to be a common feature of 
Q || magnetosheath turbulence in regimes where B rms ~ 1. a) Intensity plot of density. Only 
a segment of the simulation is shown. The presence of upstream waves is clearly evident. In 
the magnetosheath, current sheets and magnetic island can also been seen, b) A close up of Q || 
magnetosheath using LIC to show magnetic field lines colored by B. Many magnetic islands are 
observed at the shock surface all the way to the vicinity of the magnetopause. 
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FIG. 7: Plot of the total ion temperature for run 1 at £l c it — 694. Magnetic field lines are overlaid 
in white. The enhanced nragnetosheatli heating in the Q y as compared to Q± magnetosheath is 
clearly evident. Although temperature is a positive quantity, the range of color bar is only to 
improve the contrast in the figure. 
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FIG. 8: THEMIS measurements between October 2007 and October 2013 of the mean value of the 
total ion temperatures in the magnetosheath during a Parker-Spiral IMF in the Magnetosheath 
Interplanetary Medium (MIPM) reference frame. The physical dimension of each bin is 0.5 x 0.5 I?e 
and the bin number density is typically a hundred data points per bin. Each datapoint is calculated 
from the mean averaged 3-minute intervals of THEMIS data. The MSH downstream the quasi- 
parallel bow shock region is about 10-15 percent hotter. 
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FIG. 9: Formation of a hierarchy of magnetic islands due to Q || magnetosheath turbulence in a 
global fully kinetic simulation for radial IMF conditions. Only a zoomed in view of the Q y region 
is shown. The LIC of magnetic field is colored by density. 
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FIG. 10: The interaction of upstream wavefronts with the bow shock (run 5) is shown to consist 
of four stages. The waves are formed due to the relative streaming of shock reflected ions and the 
incoming solar wind. Due to the high solar wind speed, the waves get convected back towards the 
bow shock. Some of the waves steepen prior to hitting the shock (panels a-d). As the wavefronts 
go through the shock, they get amplified and compressed to shorter wavelengths (panels e-h). 
The presence of a coherent series of wavefronts inside the magnetosheath is evident in panels i-1. 
After some time, however, the wave structures get dissipated and are disrupted (panels m-p). Thus 
shortly after a radial IMF turning, coherent wavefronts should be observable in the magnetosheath, 
but at later times, only turbulence would be observed. 
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FIG. 11: In cases where the steepened wavefronts become nearly parallel to the Q y segment of the 
bow shock, they can significantly modify the bow shock and the size of the magnetosheath. The 
process stages of evolution starting fgrom impingement to eventual dissipation of the wavefronts 
indicated in Figure 10 also occur here as indicated by tracking of one wavefront labeled as WF1. 
Note the formation of a second wavefront, WF2, in panels d and h, that will repeat the process. 
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FIG. 12: Wave amplification across the shock. A zoomed in view of figure llg is shown in panel 
a. A ID plot of By along a cut in panel as is shown in panel b. Note how the upstream wave is 
amplified as it crosses the bow shock. As a reference, the red curve shows an exponential curve 
through the B y profile. 


40 



7 


0 



FIG. 13: Properties of the region in-between the shocklet and the bow shock. Close up of a shocklet 
in Figure lib is shown, a) LIC of B colored by density. The shocklet has many properties similar 
to a fast shock, showing an increase in density and magnetic field within it. Note the formation 
of magnetic islands in the region between the shocklet and the bow shock as well as along the 
surface of the magnetopause, b) LIC of ion flow colored by Vi x . Shocklet reduces the incoming 
flow speed and deflects it. An observational signature is the prediction of a turbulent solar wind 
density region upstream from the bow shock with almost no flow speed. A careful examination of 
this figure reveals presence of small scale flow vortices. 
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FIG. 14: Plot of the magnetosonic Mach number M s for run 1. The reflected ions in the fore- 
shock act to significantly lower the M s upstream of the shock. Note the high M s region in the 
magnetosheath, indicating the location of the jets seen in Figure 7. 
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FIG. 15: Generation of jets inside the magnetosheath (run 1). a) Plot of dynamical pressure in an 
area of the simulation zoomed in around the quasi-parallel bow shock. Some of the jets reach the 
magnetopause while others terminate closer to the bow shock, b) LIC of magnetic field colored by 
dynamical pressure for an area zoomed in around a strong jet marked by the bow wave in panel 
“a” that has triggered a flux transfer event at the magnetopause. 
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FIG. 16: Properties of the jets (run 1) as shown in plots of a) dynamical pressure, b) total magnetic 
field, c) total ion temperature, and d) ion density. The jets exhibit high dynamical pressure, 
enhanced magnetic field and density but cooler temperature than their surrounding plasma. The 
areas between the jets have lower density but hotter plasma in this case. 
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FIG. 17: a) Plot of Vi x for the same simulation time as in Figure 16. Sunward flow regions are 
associated with the hot regions between the jets in Figure 16c. b) The ID cut shows that the 
anti-sunward flow in the jets can have speeds comparable to the solar wind whereas the sunward 
flows can be as high as 2 Va- 
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FIG. 18: Sunward flow formation - evolution of the flow associated with the jets is shown with a 
particular focus on the flow around the bow waves (run 1). Black contours are of density. There 
is no evidence that interaction of the bow waves/jets with the magnetopause caused a strong 
sunward flow. Bow waves push the plasma and force it to go around it. This can cause sunward 
flow. Another way that the bow wave can create a sunward flow is if the jet starts to wrap akin to 
a vortex. This is evident in panel d). 
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FIG. 19: Sunward flow formation - evolution of the flow associated with the jets with a particular 
focus on the flow outside of the jets (run 1). Black contours are of density. The jets push the plasma 
downward into the lower density regions between them as evident from the blue streaklines. The 
shape and the distance between the jets change on rapid timescales as evident from tracking of 
two jets, labeled as J\ and J 2 . The flow in the region between the jets can be quite complex. 
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FIG. 20: Close up view of the two jets labeled in Figure 19 as J\ and J 2 . Note the presence of a 
region of fairly uniform magnetic field between the two jets, marked as Rl. LIC of magnetic field 
colored by ion temperature (a)-(b) and (c)-(d) V{ x at Q. c it — 736 and 756. A few magnetic field 
lines are drawn in white and the black contours are of density. The bow shock surface is clearly 
rippled. The solar wind plasma penetrates the magnetosheath almost unimpeded along the jets 
(Fig. 20c-d) whereas there is strong heating across the shock between the jets as seen in Figure 
20a-b. Note the development of sunward flow as evident from blue regions in Figure 20a-b. It is 
mainly driven by the downward flow caused by the jets (blue streaklines in Fig. 19) which push 
out the plasma between the two jets R2. 
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FIG. 21: Nonlinear formation of a large vortex (run 7). Gold curves are magnetic field lines. 
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FIG. 22: Development of the vortex in Figure 21 as seen from streaklines. Similar to Figure 19, 
there is a downward flow (pink streamlines) that pushes the plasma out sunward. In this case, the 
sunward flow gets caught in the velocity shear in the magnetosheath and wraps up into a vortex. 
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FIG. 23: a) An example of image LIC. b) An example of surface LIC. The utility of surface LIC 
is illustrated using data from a 3D fully kinetic simulation of magnetic reconnection. The 2D 
projection shows the location of the tearing “islands” which are really flux ropes as evident from 
the 3D image. The surface of the LIC is defined by isosurface of density and it is colored by electron 
speed U e . 
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